Zigzag edge modes in Z2 topological insulator: reentrance and completely flat 

spectrum 
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The spectrum and wave function of helical edge modes in Z2 topological insulator are derived on 
a square lattice using Bernevig-Hughes-Zhang (BHZ) model. The BHZ model is characterized by 
a "mass" term M{k) = A — Bk^. A topological insulator realizes when the parameters A and B 
fall on the regime, either < A/B < 4 or 4 < A/B < 8. At A/B — 4, which separates the cases 
of positive and negative (quantized) spin Hall conductivities, the edge modes show a corresponding 
change that depends on the edge geometry. In the (1, 0)-edge, the spectrum of edge mode remains 
the same against change of A/B, although the main location of the mode moves from the zone 
center for A/B < 4, to the zone boundary for A/B > 4 of the ID Brillouin zone. In the (1, 1)- 
edge geometry, the group velocity at the zone center changes sign at A/B = 4 where the spectrum 
becomes independent of the momentum, i.e. flat, over the whole ID Brillouin zone. Furthermore, 
for A/B < 1.354, the edge mode starting from the zone center vanishes in an intermediate region 
of the ID Brillouin zone, but reenters near the zone boundary, where the energy of the edge mode 
is marginally below the lowest bulk excitations. On the other hand, the behavior of reentrant mode 
in real space is indistinguishable from an ordinary edge mode. 



I. INTRODUCTION 

Insulating states of non-trivial topological order have 
attracted much attention both theoretically and experi- 
mentally. A topological insulator has a remarkable prop- 
erty of being metallic on the surface albeit insulating 
in the bulk. Recently much focus is on a specific type 
of topological insulators, [H which are said to be "Z2- 
nontrivial". Q The latter occurs as a consequence of 
interplay between a specific type of spin-orbit interac- 
tion and band structure. @ Such systems are invari- 
ant under time reversal and show Kramers degeneracy. 
From the viewpoint of an experimental realization, the 
original idea of Kane and Mele (KM) 0, Q was often 
criticized for being unrealistic, since it relies on a rel- 
atively weak spin-orbit coupling in graphene. In order 
to overcome such difficulty, Bernevig, Hugues and Zhang 
(BHZ) proposed an alternative system which is also 
Z2-nontrivial but not based on graphene. BHZ model 
was intended to describe low-energy electronic properties 
of a two-dimensional (2D) layer of HgTe/CdTe quantum 
well. Conductance measurement in a ribbon geometry Q 
showed that the system exhibits indeed a metallic surface 
state, which is also called helical edge modes. 

This paper highlights the spectrum and wave func- 
tion of such helical edge modes in the BHZ model. Re- 
specting appropriately the crystal structure of original 
HgTe/CdTe quantum well, one can safely implement it 
as a tight-binding model on a 2D square lattice. [l|, Q 
Note, however, that in contrast to KM model which can 
be represented as a purely lattice model, in BHZ an in- 
ternal spin- 1/2 degree of freedom, stemming from the 
s-type Fg and p-type Fg orbitals, resides on each site of 
the square lattice in addition to the real electronic spin. 
We also mention that in the continuum limit with van- 



ishing topological mass term, KM model has two valleys 
{K and K'), whereas BHZ has a single valley (at F). An- 
other idea which we can borrow from graphene study is 
the sensibility of edge spectrum on different types of edge 
structure, i.e., either zigzag of armchair type. [6|, |7| This 
applies also to the helical edge modes of BHZ topological 
insulator in a ribbon geometry, since the edge spectrum 
is predominantly determined by how the 2D bulk band 
structure is "projected" onto the ID edge axis. Indeed, 
the structure of BHZ helical edge modes has been ex- 
tensively studied in Ref. [1] in the (l,0)-edge geometry 
and in the tight-binding implementation. However, in 
the practical experimental setup, @ this is certainly not 
the only one which is relevant to determine the trans- 
port characteristics at the edges. Here, in this paper our 
main focus is on the other representative geometry, the 
(1, l)-edge. In the (1, l)-geometry, as a consequence of 
the specific way how " hidden" Dirac cones (or valleys) in 
the 2D spectrum is projected onto the (l,l)-axis, edge 
modes show some unexpected behaviors. 

In order to motivate further the present study, let us 
first recall the importance of edge modes in the quan- 
tum Hall state under magnetic field that exhibits a finite 
and quantized (charge) Hall conductivity a^y. In realis- 
tic samples with a boundary, quantization of Hall con- 
ductivity is attributed to dissipationless transport due to 
a gapless chiral edge mode. In contrast to charge Hall 
effect, a finite spin Hall current does not need breaking 
of the time reversal symmetry. In the quantized spin Hall 
(QSH) effect, the Chern number in the bulk takes an in- 
tegral value, and correspondingly there appears integral 
pairs of gapless edge modes. On the other hand, the Z2 
topological insulator is characterized by an odd number 
of gapless modes per edge that are robust against weak 
perturbations preserving the time-reversal symmetry. 



The existence of such gapless edge states is gener- 
ally guaranteed by a general theorem under the name of 
bulk/edge correspondence. [1,0] The BHZ model has a 
convenient feature that the location of the minimum en- 
ergy gap can be controlled by changing the parameters 
in the model. In particular, the sign of spin Hall conduc- 
tivity changes discontinuously as the mass parameter of 
the model is varied. Hence, the corresponding change of 
edge spectrum should provide useful information on the 
bulk/edge correspondence. Furthermore, understanding 
of the nature of edge modes under specific edge geome- 
tries should serve as possible application of topologically 
protected phenomena in nano-architectonics. We take 
the representative cases of the (1,0)- and (l,l)-edges, 
which we call also the straight and zigzag edges, respec- 
tively. 

This paper is organized as follows: In Sec. II, we 
clarify our motivations to study the lattice version of a 
Zi topological insulator (BHZ model), implemented as 
a nearest-neighbor (NN) tight-binding model. Explicit 
form of the BHZ lattice hamiltonian is introduced in 
Sec. III. It is demonstrated that by considering a lat- 
tice model, one naturally takes into account hidden Dirac 
cones, the latter lacking in the analyses based on the con- 
tinuum Dirac model. In Sees. IV and V, we study the 
detailed structure of gapless edge modes under two dif- 
ferent edge geometries: straight and zigzag edges. Sec. 
V is the highlight of this paper, demonstrating that the 
zigzag edge modes of BHZ lattice model show unique 
features. We first point out that at A = AB, a pair of 
completely flat spectrum appears at i? = 0; besides the 
edge wave function can be trivially solved. We then show, 
in a half-empirical way, that this analytic exact solution 
at A = 4i? can be generalized to the case of an arbi- 
trary A/_B e [0,8] (this idea is schematically represented 
in FIG. H]). Using the exact solution thus constructed, 
we highlight the nature of reentrant edge modes, another 
unique feature of the edge modes in the (1, l)-edge geom- 
etry. The reentrant edge modes possess two contrasting 
characters in real and momentum spaces: though well 
distinguished in real space, they live in an extremely 
small energy scale in the spectrum. Sec. VI is devoted to 
conclusions. The gapless edge modes of BHZ topological 
insulator are also treated in the framework of continuum 
Dirac model in Appendix A. 



II. STATEMENT OF THE PROBLEM 

It has been well recognized that the quantized spin Hall 
conductivity is determined by wave functions of Bloch 
electrons over the entire Brillouin zone. On the other 
hand, only the electronic states near the Fermi level are 
relevant to the change of the Hall conductivity when the 
Fermi level is shifted. Simplified effective models are use- 
ful for the latter case since various theoretical techniques 
can be employed in the low-energy range. In this pa- 
per, we work mainly with the lattice version of the BHZ 
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FIG. 1: Recipe (conceptual) for constructing the exact edge 
wave function in the zigzag edge geometry — a half-empirical 
way. 



model, and make some comments in the low-momentum 
limit. 



A. Continuum vs. lattice theories 

Why do we have to go back to a lattice model? Firstly, 
because we need to recover the correct absolute value of 
spin Hall conductance a^y. The latter is defined as the 
difference of Hall conductance for up and down (pseudo) 
spins multiplied by ~h/{2e): 



-Ye^< 



(1) 



In quantized spin Hall (QSH) systems, the spin Hall con- 
ductance is quantized to be an integer in units of e/(27r). 
[26j This is completely in parallel with the quantization 
of charge Hall conductance in units of /h in quantized 
Hall system. In both cases, such integers are topological 
invariants and protected against weak perturbations. 

On the other hand, if one calculates, using Kubo for- 
mula, the contribution of a single Dirac fermion, e.g., of 
the continuous Dirac model at the F-point [4] to spin Hall 
conductance, then this gives half of the value expected 
from the topological quantization. [fSj - TlGf In order to 
be consistent, it is naturally assumed that there must 
be even number of Dirac cones. jl7t However, the low- 
energy effective theory with which we are starting con- 
tains obviously a single Dirac cone. 0,111 As we will see 
explicitly, such trivial discrepancy is naturally resolved 
by considering a lattice version of the BHZ model. 

Another aspect motivating us to employ the lattice ver- 
sion of BHZ model is the fact the idea of an edge states 
is a real space concept, and we need a priori to go back 
to real space to give an unambiguous definition to it. In 
this paper, we highlight the detailed structure of gap- 
less edge modes under a specific edge geometry. Clearly, 
introduction of the latter needs a precise description in 
real space. Recall also here that edge modes of graphene 



3 



nano-ribbon exhibit contrasting behaviors in zigzag and 
armchair edge geometries: e.g., the system becomes ei- 
ther metalhc or semi-conducting in the armchair geome- 
try, depending on Nr mod 3, with iV^ being the number 
of rows, whereas a completely flat edge mode appears 
in the zigzag geometry, p, 01 Where does the difference 
comes from? In momentum space, the question is how 
the bulk Dirac cone structure look like when viewed from 
the edge. Note that in the zigzag geometry the flat edge 
mode connects the two Dirac points: K and K' , whereas 
in the armchair geometry these two points are projected 
onto the same point at the edge. In a topological insu- 
lator, this bulk to edge projection is even a more subtle 
issue, since not all the Dirac cones are explicit (see Ta- 
ble I). In a sense, zigzag edge in the square lattice BHZ 
model (see FIG. [S]) plays the following double role: it re- 
sembles the zigzag edge in graphene at A = 4i?, whereas 
it may rather resemble the armchair edge at A = and 
at A = 8B. 



B. Continuous Dirac model and its boundary 
conditions 

Spin Hall conductance is a topological quantity deter- 
mined by the global structure of entire 2D Brillouin zone. 
The helical edge modes, encoding the same information, 
is, therefore, a priori not derived from a local description 
in the 2D Brillouin zone. One exception to such a gen- 
eral idea is the study of Ref. [l^l (see also Appendix) , in 
which gapless edge modes are derived from the contin- 
uous model in a strip geometry by simply applying the 
condition that all the pseudo spin components of wave 
function vanish at the boundary. This implies that 
information about the helical edge modes is actually en- 
coded in the original (single) Dirac cone. This seems 
to be rather surprising, if one recalls that in the case 
of KM model, the distinction between trivial and non- 
trivial phases is made by a relative sign of the mass gap 
at K- and K'- points, which are, of course, macroscop- 
ically separated in momentum space. Here, in the BHZ 
model, the same distinction is made by relative sign be- 
tween the mass (fc°-) term and the k^- term added to the 
Dirac Hamiltonian. 

Motivated by this observation, we investigate the 
structure of helical edge modes under different bound- 
ary conditions for the periodic BHZ model, implemented 
as a square lattice and nearest-neighbor (NN) tight- 
binding model model. In parallel with the arm chair and 
zigzag edges for graphene, we consider (a) usually consid- 
ered (1,0)- (straight) boundary as well as (b) (1, 1)- 
(zigzag) boundary for the tight-binding BHZ model. 



C. Explicit vs. hidden Dirac cones 

The idea of focusing on Dirac fermions in the descrip- 
tion of quantized Hall effect has appeared in the context 



of transitions between different plateaus. For describing 
the transitions, half-integer quantization is not a draw- 
back, since the difference of Hall conductance before and 
after the gap closing is quantized to be an integer in 
units oie^/h: 1/2 — (—1/2) = 1 or vice versa. A discrete 
jump in the Hall conductance across the transition is in- 
deed consistent with counting based on the emergence of 
massless Dirac fermions at the transition [11, 21J. 

The absolute value of Hall conductance is, on the other 
hand, a winding number, and determined by the vor- 
tices [l3|. Here, each vortex gives, in contrast to a Dirac 
cone, an integral contribution to the Hall conductance 
(in units of e'^/h). An interesting question is whether 
the total Hall conductance, often expressed as a topo- 
logical invariant [lo| . can be also written as a sum of 
contributions from emergent Dirac electrons in the spec- 
trum. Our empirical answer is yes, [17] indicating that the 
number of Dirac electrons emergent in the spectrum is al- 
ways even, reminiscent of the no-go theorem of Nielsen 
and Ninomiya in 3-1-1 dimensions [l^. It should be noted 
that here not only explicit Dirac electrons (gapless for a 
given set of parameters) but also hidden Dirac electrons 
(massive for that set of parameters) must be taken into 
account. Such massive Dirac electrons are called "spec- 
tators" in Ref. [l^, in the sense that they are inactive 
for the transition. Spectators are indispensable to ensure 
the correct integer quantization of the Hall conductance. 
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FIG. 2: Bulk energy spectrum: E = E(kx, ky) (vertical axis) 
of a BHZ lattice model, here, implemented as a NN tight- 
binding model on a square lattice: cf. Eqs. @. The spectrum 
is shown over the entire Brillouin zone: — vr/a < k^, ky < vr/a 
(horizontal plane). Parameters are chosen such that A = B = 
1 and A = 2, i.e., the system is in the topologically non-trivial 
phase: < A/B < 8. The spectrum shows a typical wine- 
bottle structure around the F-point. 



III. BHZ MODELS 
A. BHZ model in the long-wave-length limit 

Let us first consider the BHZ model in the long-wave- 
length limit. The low-energy effective Hamiltonian, de- 
scribing the vicinity of gap closing at F = (0,0), is the 
minimal model to capture the physics of a .Z2-topological 
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TABLE I: Nature of four Dirac cones in the BIfZ lattice model. The four Dirac cones appear at different values of the tuning 
parameter A, and at different points of the BZ: F, X, X' and M. Away from the gap closing, such Dirac electrons acquires 
a mass gap. The sign of such mass gap, together with the chirality x, determines their contribution to a'^y — ±e^//i. In the 
table, only their sign is shown. The symmetry of the valence orbital is also shown in the parentheses, which is either, s (inverted 
gap) or p (normal gap), corresponding, respectively, to the parity eigenvalue: 5a = +1 or (5p = — 1. The latter is related to Z2 
index u as (—1)" = IIdp ^dp- [ISI We also assume B > with no loss of generality. 
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insulator. This effective Hamiltonian is also contrasting 
to the prototypical KM model, in that the former de- 
scribes only a single Dirac cone. The distinction between 
the Z2- trivial and non-trivial phases is, therefore, made 
by adding a /c^-tcrm to the usual Dirac Hamiltonian. The 
explicit form of BHZ Hamiltonian is implemented by the 
following 4x4 matrix: 



H{k) 



hik) 
h*i~k) 



(2) 



where k — [k^, ky) is a 2D crystal momentum, here mea- 
sured from the F-point. The lower-right block h*{—k) is 
deduced from h{k) by imposing time reversal symmetry. 

The bulk energy spectrum: E = E{k) is then given by 
solving the eigenvalue equation for the upper-left block 
hik) of the 4x4 BHZ Hamiltonian, i.e., 



hik)'iPik) = EiPik). 



(3) 



In order to represent hik) in a compact form, we in- 
troduce a d- vector, d — idx,dy,dz), each component 
of which is either an even or an odd function of fc: 



''X,y,z 



(k), whose parity is determined by symme- 



try considerations. [J] As far as the low-energy universal 
properties in the vicinity of T-point is concerned, we need 
to keep only the lowest order terms of dx,y,zik), and in 
this long-wave-length limit, they read explicitly. 



dxik) — Akx, dyik) — Aky 
4(fc) = A-Bikl + kl). 



(4) 



Other parameters which appear in Ref. i.e., C and 
D are set to be zero, which, however, does not lose 
any generality. The bulk energy spectrum is thus de- 
termined by diagonalizing the following "spin Hamilto- 
nian", hik) = d(fc) • cr, where cr = {gx, (Jy, (Jz), are Pauli 
matrices. Using the standard representation for tr, hik) 



reads explicitly as. 



hik) = 



dx 



(5) 



Each row and column of Eq. ([3]) represent an "orbital 
spin" associated with the s-type Tg and the p-type Tg 
orbitals of the original 3D band structure of HgTe and 
CdTe. [1^ Then, by choosing the "spin quantization 
axis" in the direction of d-vector, one can immediately 
diagonalize the Hamiltonian /i(fc), i.e., 

hik)\dik)±) ^ ±dik)\dik)±) , (6) 

where the eigenvalue Eik) = ±(i(fc) given by. 



dl 



dl 



dl. 



(7) 



|fcp+B2|fc|4 

> A^I2B = 



, (8) 
Aq, the 



dik) = |d(fc)| = 

This implies, 

Eik)"^ = A2 + {A^ - 2BA) 

where IfcP = kl + kf,. When A 
dispersion relation ^ represents a wine-bottle structure 
(FIG. 2), i.e., i?(|fe|) shows a minimum at a finite value 
of |fc|. At the critical value Aq = A'^/2B, the density of 
states shows van Hove singularity. 



B. BHZ model on square lattice and Dirac-cone 
interpretation 

Lattice version of the BHZ model is implemented as a 
tight-binding Hamiltonian. To construct such a Hamil- 
tonian explicitly, we replace linear and quadratic depen- 
dences in hik) on kx and ky as in Eqs. (O, by a function 
which has the right periodicity of the square lattice. This 
can be implemented as. 



dxik) 
dyik) 



dzik) A 



2B 



j- sinikxa), 
^ sinikya), 
[2 — cos(fc2;a) — cos(A:j^a)] , (9) 
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FIG. 3: ai^J in units of e/2Tv, obtained by numerically eval- 
uating the fc-integral in Eq. pip, plotted as a function of 
A/B. 



where a is the lattice constant. Eqs. ([9]) corresponds to 
regularizing the effective Dirac model on a square lat- 
tice with only nearest-neighbor (NN) hopping. In this 
setup, i.e., Eqs. (P together with Eqs. (P, the 

lattice version of BHZ model acquires four gap closing 
points shown in TABLE I, if one allows the original mass 
parameter A to vary beyond the vicinity of A = 0. The 
new gap closing occurs at different points in the Brillouin 
zone from the original Dirac cone (F-point), namely at 
Xi = (7r/a,0), X2 = (0,7r/a) and M = (n/a^Tr/a). The 
gap closing at M occurs at A = 85, whereas the gap clos- 
ing at Xi and X2 occurs simultaneously when A = AB. 

Each time a gap closing occurs, one can re-expand the 
lattice model with respect to small deviations of k mea- 
sured from the gap closing. The new effective model 
in the vicinity of such hidden gap closing falls on the 
same Dirac form as the original one at the F-point, up 
to the A:^-term. In order to quantify the emergence of 
such hidden Dirac cones, one still needs the following 
two parameters: (i) the mass gap A (especially, its sign), 
and (ii) the chirality x- The latter is associated with 
the homotopy in the mapping: k d{k). Note that in 
the gap closing at Xi and X2 at A = AB, the role of 
k± = ± iky is interchanged compared with the origi- 
nal Dirac cone at the F-point. The former (latter) corre- 
sponds to X = ~1 (x = +!)• The missing Dirac partner, 
in the sense of Ref. [l^ , is found in this way. Once the 
explicit form of effective Dirac Hamiltonian in the con- 
tinuum limit is given, one can determine its contribution 
to criy"''. In systems with TRS, i.e., of the form ([2|), the 
contribution from h{k) to criy — a^y + ajiy cancels with 
that of h*{~k). On the other hand, their contribution 

to crit/ = —(e/2'K)(a'ly — a^y), remains finite, takes a 
half-integral value, ±1/2 in units of e/27r (2e^//i in con- 
ductance). Contribution to aiy from a Dirac point with 
a mass gap A and chirality x is, 

4;)=Xsign(A) i ^, (10) 



provided that the Fermi energy is in the gap. This can 
be verified explicitly by applying the Kubo formula to 
the continuum model. 

As mentioned earlier, such counting based on the con- 
tinuum Dirac model, is known to describe correctly a 

(c) 

discrete jump of axy in the quantum Hall case. Here, 
we apply the same logic to QSH case. Imagine that one 

(s) 

observes the evolution of axy , starting with the trivial 

(s) 

insulator phase, where axy = 0, and varying the mass 
parameter A. The Fermi energy is always kept in the 
gap unless there appears a Dirac cone. Each time such 

(s) 

a gap closing occurs, axy shows a discrete change, which 
can be attributed to the above Dirac fermion argument. 
In the present model, one can verify explicitly that this 
is indeed the case. 

If one evaluates the spin Hall conductance from TKNN 
formula, ai^y allows for the following representation, 
in terms of Berry curvature integrated over the entire 
Brillouin zone: 

e P (Pk dd dd d 

where d = |<i| and d = d{k) is given, e.g., by Eq. 
For such an explicit choice of d{k), Eq. (fTTj) is evaluated 
numerically, and plotted as a function of A/i? in Fig. |31 
When d{k) is given by Eq. the plotted curve (the 
solid curve shown in blue in Fig. [3l which looks practi- 
cally like steps) is comparable with the column X]_dp 

of TABLE I. Note that the absolute value of axy is sus- 
ceptible of the concrete implementation of d{k) over the 
entire Brillouin zone, whereas its parity (whether it is 
even or odd) in units of e/(27r) remains the same. As 
well known, the latter determines system's Z2-property. 

We have seen that a'^y takes a finite value ±e/27r when 
< A/B < 8, i.e., the system is in the topological (in- 
verted gap) phase. This is also consistent with the gap- 
less edge picture, in which the spin Hall conductance of 
twice the unit of quantum conductance /h is attributed 
to two channels of edge modes, which form a pair of 
Kramers partners. The apparent half-integer quantiza- 
tion at the F-point, in the sense of Eq. ([TU]). is compen- 
sated by the contribution from missing Dirac partner(s), 
and as a result, is indeed shifted by one-half, replaced by 
an expected integral quantization. 



IV. STRAIGHT EDGE GEOMETRY 

Let us first review the behavior of gapless edge modes 
in the straight edge geometry, the latter commensurate 
with the square lattice, and can be chosen either nor- 
mal to the (1,0)- or to the (0, l)-direction (as in FIG. 
U). Introducing an edge leads to breaking of the trans- 
lational invariance in the direction perpendicular to the 
edge, inducing a coupling between Dirac cones. 
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FIG. 4: Straight edge geometry. The two boundaries of the 
strip (=edges) are, here, chosen to be perpendicular to the 
(0, l)-direction. In this figure, the number of rows in the strip 
is iV^ = 5. 



the edge. Eq. (IT3|) thus rewrites, 

(15) 

where Z?(fca;)'s are diagonal (on-row) components, which 
read explicitly, 

D{kx) = AcTa; sin fca; + {A - 2B (2 - cos fca;)} (Tz 

= AaxSinkj: + Q{kx)(Tz- (16) 

Tj^i^j represents a hopping amplitude in the y-direction, 
i.e., between neighboring rows. Inside the strip, i.e., for 
J — 1, • • • ,Nr, these amplitudes take the same value as 
in the bulk, given in Eqs. ([T4| . i.e.. 



r 



.1+1, J 



= Ty = -i- 



A 



(17) 



A. Effective one-dimensional model 

In the straight edge geometry shown in FIG. |4|), elec- 
trons are confined inside a strip between the rows a.ty — a 
and y = Nj-a. The translational invariance along the x- 
axis is still maintained, allowing for constructing a ID 
Bloch state with a crystal momentum k^'. 



J), 



(12) 



where k^ is measured in units of 1/a with a being the 
lattice constant. |/, J) = c| j|0) is a one-body electronic 

state localized on site (/, J), and c| j is an operator cre- 
ating such an electron. It is also convenient to introduce 
c|. J, and express \kx, J) as \kx, J) — cl j|0). Naturally, 
the two creation operators are related by Fourier trans- 
formation similarly to Eq. (|T2)) . i.e., c|. j = J2i e^^^c\ j. 

In order to introduce the edges, it is convenient to 
rewrite the BHZ tight-binding Hamiltonian in terms of 
the hopping between neighboring rows. Let us first con- 
sider the BHZ tight-binding Hamiltonian in real space: 

iJ = ^ [(A - AB)<j4 jCj,j 
i,j 

+ |r:Ec|_^i_jC/,J -f Fycl j^-^C/^J + /l.C.j ,(13) 

where 2x2 hopping matrices, F^; and Fj,, are given ex- 
plicitly as. 



.A 



A 



= -i—ax + BcTz, Ty = -i—<7y + BcTz 



(14) 



In order to rewrite it in terms of the Bloch state, Eq. 
(I12p . or equivalently, in terms of the corresponding cre- 
ation and annihiration operators, cj. j (ck^j), we per- 
form Fourier transformation in the (x-) direction along 



In the tight-binding implementation, a strip geometry 
can be introduced by switching off all the hopping ampli- 
tudes connecting sites on the edge of the sample to the 
exterior of the sample. In our straight edge geometry, 
such outermost rows are located at J = 1 and J = Nr. 
We turn off all the hopping amplitudes from J = 1 to 
J = 0, and the ones from J ~ Nr to J ^ Nr + 1, i.e.. 



0. 



(18) 



This boundary condition, (i) breaks the translational in- 
variance in the y-direction, and (ii) restricts the Hamil- 
tonian matrix into Nr x N^. blocks. 



B. Spectrum and wave function 

Let us construct the eigenvector of the straight edge 
Hamiltonian, Eq. (llSp . with an eigenenergv E. Since Eq. 
(IT5|) is already diagonal w.r.t. kx, we diagonalize Eq. 
for a given k^, to find the energy spectrum E — E{kx). 
The corresponding eigenvector is thus specified by E and 
kx, and takes generally the following form: 



\E,kx} 



XT J) J 



(19) 



where ipji-l^j kx) is a 2 x 2 spinor specifying the amplitude 
and the pseudo spin state of eigenvector on row j. One 
might rather regard, 



-02 
V'3 



(20) 



as the wave function of the corresponding eigenstate. The 
eigenvalue equation. 



H\E,k.x)^E{kx)\E,kx) 



(21) 
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kx/;r kal K kx / 7i 

FIG. 5: Energy spectrum (numerical) in the straight edge geometry for different values of A (j4 = B = 1). The number of 
rows Nr is, here, chosen to be Nr — 100. The dotted curve is a reference, showing the exact edge spectrum given in Eq. (|33|l . 
Starting with the left panel A = B (spectrum shown in red), A — 4B (center panel, spectrum in green), and A — 5B (right, 
blue). 
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FIG. 6: \pi\ and \p2\ plotted as a function A/B in the limit 
0. 



can be rewritten, in terms of the ipj{E, k^Ys, in the form 
of a recursive equation: 



(22) 



All the information on the spectrum and the wave func- 
tion of both the extended bulk states and the localized 
edge states is encoded in Eq. ([^ and the boundary 
condition which we will specify later. Since the recursive 
relation, Eq. ([22l) . is linear, its eigenmodes take the form 
of a geometric series: 



(23) 



where p is a solution of the characteristic equation which 
we will derive later, li \p\ < 1, Eq. (l23l) may represent 
an edge mode. Since the recursive relation, Eq. ([22| . 
is of second order, its characteristic equation becomes a 
quadratic equation, giving two solutions for p, say, p — 
Pi,2- On the other hand, our recursive equation has also 
a 2 X 2 matrix form, we first have to solve a (reduced) 
eigenvalue equation for ^E'o, assuming that p is given. The 
reduced eigen value equation for reads. 



D{k^)+pTy 



irt 

p \ 



ipo = Eipo. 



(24) 



Using Eqs. (ITi)) this can be also rewritten as, |2S 



D{k^) + i 



A 



- pj ay + B \^-+ pj az -00 = EtPo, 

^ (25) 
where D{kx) is given in Eq. This is a 2 x 2 eigen- 

value equation, and there are generally two solutions for 
E and two corresponding eigenvectors for a given p. Re- 
call here that for < A < 4, our numerical data (FIG. 
El) show that the edge spectrum behaves as i? — in the 
limit of kx — > 0. Hereafter, we will focus only on such 
edge solutions. Since in the same limit, D{kx) — >■ ^{0)az, 
Eq. ([25]) reduces to. 



A 
2" 



B 



P 



00 = 0. (26) 



Note that we have multiplied both sides of Eq. (1^ by 
az- It is clear from this expression that \E'o can be chosen 
to be an eigenstate of i.e., ^'q — |a;±), where 



1 

71 



1 

71 



1 

-1 



(27) 



If one denotes the eigenvalue of Ux by s, as Ux'^q = 
±5*0 = s^fQ, then s = 1 corresponds to ^'o = \x+)^ 
and s = — 1 to = \x—). Namely, s specifies the eigen- 
spinors given in Eqs. (j27p . Based on these eigenspinors, 
one can construct the total wave function. Of course, one 
still needs to do determine the allowed values of p. For 
an eigenstate specified by s, p must satisfy. 



17(0) 

For s = 1, 



A 
2" 



B 



P 



0. 



(28) 



-O(O) ± Vf7(0)2 -4^2 _ 



A + 2B 



= Pi,2(0). (29) 



In Eq. p8| , if p is a solution of this quadratic equation for 
s = 1, then 1/ p satisfies the same equation for s = — 1. [l| 
Thus the general solution becomes a linear combination 
of the following four basic solutions: 



^3 



[c+iPiW 
C-iPi(O)- 



■^+2P2W] 
C-2P2(0)" 



(30) 
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Of course, at this point, this is just a solution at only 
one single fc-point, = 0. However, a solution in the 
form of Eq. (l30|), with pi,2(0) given in Eq. ([29|), can 
be easily generalized to satisfy Eq. (1^^ for an arbitrary, 
finite kx, after a simple replacement of parameters. As for 
the eigenmode of the form of Eq. ([23l) , one has to solve 
a reduced eigenvalue equation for ipQ, and determine p 
such that Eq. (I^SI) is satisfied. However, since D{kx) has 
the structure given in Eq. (ITG]) . the eigenmodes for ijjQ 
given as Eq. (j27p remain to be valid for an arbitrary fc^. 
This might become clearer, if one decomposes Eq. (P5|) 
into the following set of equations: 



Aa^sinfc^^'o = -B^'o 



.P 



*n = 0, 



(31) 



(32) 



i.e., Eq. (I?5]) is recovered by adding both sides of Eqs. 
([5T|) and ([5^ . The first equation gives the energy dis- 
persion, E = Es{kx), if "00 is chosen to be an eigenstate 
of ax, i.e., cra;V'o = ±"00 = stAo, where 



Es{k.. 



sA sin k„ 



iAsin fc^. 



(33) 



Note that this is an exact edge spectrum valid over the 
entire range of k^ , as far as the edge solution is possible 
(see FIG. [S]). On the other hand, Eq. (15^ . analogous to 
Eq. (1^51) . justifies the previous conjecture: -00 = |2;±). 
While, in the characteristic equation for p, one has to 
make the simple replacement : fi(0) — )• ^{k^), i.e., Eqs. 
([5^ and ([^5]) are identical up to this replacement. For 
s = 1, the solution for p reads, 



P 



-n{k.,)±^n{kxY + A^-AB^ 
A + 2B 



PiAkx). (34) 



Correspondingly, a general solution for ipj can be con- 
structed as, 

0j = [c+ipi{kj;y + c+2P2{kxy] 

+ [c-ipiik^y^ +c^2P2{kxy^\x-), (35) 

where the coefiicients c±i_2 should be chosen to satisfy 
the boundary conditions. Eq. (PSj) is smoothly connected 
to Eq. ^ in the limit: k^ 0. 

C. Illustration of edge spectrum 

Three panels of FIG.[5]show the energy spectrum (edge 
-I- bulk) for different values oi A/B. As for the edge part 
of the spectrum, only a part of Eq. p3p is realized. In 
order to determine which part of the spectrum in Eq. 
([33)1 is indeed activated, we discuss below the case of 
semi-infinite geometry in some detail. 

FIG. [S] also demonstrates one of another specific fea- 
ture of straight edge mode that the main location of the 



mode moves from the zone center for A < 4B, to the 
zone boundary for A > 4B. Thus, the group velocity in- 
tersecting with the Fermi level reverses its sign, reflecting 
the sign change of a^y in the bulk. This can be regarded 
as the concrete expression of bulk/edge correspondence 
in the present case. [1, |^ 

What kind of a boundary condition should we apply in 
Eq. (1351) ? Suppose that here our system is semi-infinite, 
for simplicity, extended from j — 1 to j oo. Such a 
boundary condition can be applied, by formally requiring 
that the wave function ((55)) vanishes at j = 0, i.e.. 



00 



(36) 



This means that the coefficients c±i.2 in Eq. (|35p should 
be chosen to satisfy. 



c+i + c+2 = 0, 



c_i + C_2 



= 0. 



(37) 



This turns out to be rather an important requirement for 
determining the range of validity of the solution given in 
Eq. (j35p . since the wave function ipj must be normaliz- 
ablc. In Eq. ([35]), only the eigenmodes of the form of 
Eq. ([23)) with |p| < 1 should be kept in the solution (to 
be precise, both \pi\ and \p2\ must be smaller than 1). 
In a strip geometry, another solution, consisting of both 
\p\ > 1, describes the edge mode localized at the opposite 
end of the system. 

In FIG. [HI \pi I and \p2 \ are plotted as a function of A/B 
in the limit k^ 0. When < A/B < 4, both \pi\ and 
\p2\ are smaller than 1, namely both |l/pi| and |l/p2| are 
larger than 1 . This means that only the first two terms 
of Eq. ([55)1 . both corresponding to (s = 1), should 
be kept in the solution, i.e.. 



c+i = -c+2 7^ 0, 



C-2 



0. 



(38) 



Outside this region, either \pi\ or \p2\ is larger than 1. 
When IpI > 1, since this implies automatically \l/p\ < 1, 
the eigenmode corresponding to the latter is still compat- 
ible with the boundary condition at j ^ oo. However, 
because of the boundary condition at j — 0, i.e., Eq. 
([37)) . when |pi| > 1 and \p2\ < 1, or vice versa, the only 
possible choice for the coefficients c±i^2 is 



c+i = c+2 = c_i = C_2 



= 0. 



(39) 



Namely, a solution of the type of Eq. (|35p. or an edge 
mode crossing at fc^; = is inexistent. This is consistent 
with the fact that an edge mode crossing at fc^; = ex- 
ists only in the region, A/B £ [0,4] in the straight edge 
geometry, (when A/B G [4,8], the edge modes cross at 
kx = TT, i.e., at the zone boundary, which is also time- 
reversal symmetric.) 

Coming back to the regime in which edge modes are 
existent, i.e., A/B € [0,4], one can clearly see in FIG. 
[S] that there are two different behaviors — a flat region 
where |pi| and \p2\ are degenerate, and the remaining 
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FIG. 7: Upper panel: \pi\ and \p2\ plotted as a function of 
fcj; at A = 2 (A = B = 1). Reference lines are at k^/Tr — 
±0.5, and at \p\ = 1 and \p\ = l/VS — as for the latter, cf. 
Eq. (|40|) . \p\ — 1 corresponds indeed to the point at which 
the edge modes merge with the bulk spectrum (central panel, 
Nr — 50). Lower panel: an enlarged image of the central 
panel, plotted also for a system of larger size: A^,- — 100. 



part with two branches. This is due to the fact that 
the two solutions for p could be either both real, or a 
pair of complex numbers conjugate to each other. In 
the latter case, the two solutions have the same absolute 
value, \pi\ — IP2I, whereas in the present case, one can 
verify that this degenerate value is independent of A, i.e.. 



iPil 



\P2\ 



A-2B 



(40) 



j A + 2B 

This explains the existence of a flat region in FIG. IH) 



From Eq. (l29)) we see that the square root becomes pure 
imaginary for all k provided A_ < A < A_|_, where [29| 



A± 
B 



2 1±A 1- 



452 



(41) 



At IpI = 1, the edge solution is expected to merge with 
the bulk spectrum. This happens, when 



cos kx = 1 . 

2B 



(42) 



Such a behavior becomes clearer by plotting |p|'s as a 
function of /c^.. FIG. [7] illustrates this feature at A = 2 
for A — B ~ 1. One can indeed see that the edge spectra 
merge with the bulk at k^ = km, satisfying Eq. (|42l) . 
The latter reduces, at this value of A, to cosfc^ = 1 — 
A/{2B) = 0, i.e., k„, = ±n/2. 

It is also instructive to investigate the nature of edge 
modes in real space, i.e., the wave function, and compare 
it with the general solution psp . In numerical experi- 
ments, one has to diagonalize the 2Nr x 2Nr Hamiltonian 
matrix, equivalent to Eq. (I15p . An eigen wave function 
is, therefore, obtained as a 2A'^r-component vector; here, 
in the straight edge geometry, the latter can be chosen to 
be real. The edge wave function is easily identified if it 
exists, e.g., by choosing the lowest-energy eigenmode vJ^q 
in the upper band. By investigating the structure of such 
an edge wave function, one can explicitly verify that the 
eigenmodes are spanned by two eigenspinors given in Eq. 
(|27p . In repeating such numerical experiments for differ- 
ent k^ and A/B, one can naturally distinguish an edge 
state from a bulk state by focusing on the spatial distri- 
bution of the wave function. Here, what deserves much 
attention is that one can recognize a one-to-one corre- 
spondence between localizability of the wave function ^'o 
and its spinor structure. 



V. ZIGZAG EDGE GEOMETRY 

Let us turn to the case of a different edge geometry, the 
zigzag edge geometry, shown schematically in FIG.[S] As 
mentioned earlier, the zigzag edge geometry considered 
here is, in a sense, analogous to a more popular edge ge- 
ometry of graphenc ribbon, named in the same way, but 
defined on a hexagonal lattice. Here, on a square lattice, 
a zigzag edge is introduced, either normal to (1, 1)- or 
(1, — l)-direction (as in FIG. [H]). Electrons in the zigzag 
edge geometry are, therefore, confined to a strip diag- 
onal in the cartesian coordinates, say, along the (1,1)- 
direction as in FIG.[8l 

Intuitively, say, because the zigzag surface is literally, 
"rough", one expects that this edge geometry might have 
a stronger tendency to trap electrons in the vicinity of 
the boundary. We show below, on one hand, that this 
intuition from the macroscopic world is still valid in the 
microscopic quantum mechanical world. But just as a 
result of this stronger tendency to keep the electrons in its 
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FIG. 8: Zigzag or (1, l)-edge geometry (here, chosen to be 
normal to the (1, — l)-axis), defined in terms of the original 
square lattice, on which the new x- and y-axes are superposed 
in blue (online). Numbers on these axes are redefined indices 
/ and J, i.e., (i, y) = {la/ \/2, Ja/ \/2). The number of rows 
Nr is here chosen to be, Nr ~ 6. When A'^r is even (odd), the 
two edges are inversion asymmetric (symmetric) w.r.t. the 
center of strip. 



vicinity, on the other hand, the zigzag edge shows various 
curious phenomena, e.g., completely flat edge modes, and 
the reentrance of edge modes in fc-space, etc. 

Clearly, the translational invariance along the (1, 1)- 
axis is maintained, on which i-axis is introduced, to- 
gether with the conserved momentum fc = fcj in this di- 
rection. Accordingly, y-axis is chosen to be in the (1,-1)- 
direction. It may be also useful to redefine the indices / 
and J such that the lattice points in the original square 
lattice are located at (i,y) = (/a/\/2, Ja/\/2), where 
I is an even (odd) integer for J: even (odd). In the 
zigzag edge geometry, the spectrum E = E{k) and the 
wave function are determined by the following recur- 
sive equation: 



(A - AB)a,^j + T4,,+i + T^^j.i = 



(43) 



for , analogous to Eq. (|22|) in the straight edge geom- 
etry. In order to derive Eq. (j43p . we first rewrote the 
tight-binding Hamiltonian (15) in the new labeling, and 
then considered a Bloch state along the ai-axis, analogous 
to Eq. (jl2p . but with a crystal momentum k conjugate 
to i, i.e., k = kx. As was the case in Eq. ((22|) . F de- 
scribes, here, in Eq. (j43p the hopping between adjacent 
rows, and reads explicitly as. 



r - 



-ik ^ 



■ ^ ik 

i—e a. 
2 



y 



2Ba^ cos k. 



(44) 



Note that here the crystal momentum k is measured in 
units of l/(-\/2a) so that the zone boundary is always 
given by fc = TT. (soj To find an edge solution, we first 
express the solution of Eq. in the form of a geo- 

metric series, written formally in the same as Eq. p3|. 
Recall that p is (generally) a (complex) number of, for an 
edge mode, amplitude smaller than unity (with the un- 
derstanding that the edge mode is localized in the vicin- 



ity of J = 0). ipQ is a two component eigenvector of the 
following reduced eigenvalue equation: 



(A - AB)a, +pr+ -T^ 
P 



ijo = Eijo- (45) 



From the analogy to the straight edge case, one may ex- 
press r as 



A 

r = i'-^Ck{(Tx - O-y) + -^Sk{(Ta: + ay) + 2BckCrz, (46) 



.A 

2"' 



and rewrite Eq. (j45l) into the following explicit form: 



(A - AB)a, 



P 



A 



Ski^x +(Ty) + 2Bcfe(Tj 



1^ ( ^ 



ipo = EipQ. 



(47) 



where we Ck and Sk are short-hand notations for, respec- 
tively, cosfc and sinfc. Comparing this form with the 
straight edge case, one can see that here one cannot use 
the same recipe for solving the problem, i.e., solving the 
problem at, say, k = and then extrapolate its solution 
to general k. It seems not impossible to proceed in that 
direction and solve the problem analytically, but here, we 
choose to take another route, which is much simpler, to 
find still an exact solution, but in a half-empirical way. 



A. Completely flat edge mode at A = iB 

Some concrete examples of such energy spectrum are 
shown FIG. 13 A pair of gapless edges modes always 
appear iff < A/B < 8. In contrast to the straight 
edge case, however, they appear always in the vicinity of 
k = 0, and intersects at A; = 0. 

One of the last panels of FIG. [3] shows a unique feature 
of edge modes in the zigzag edge geometry. At A = 4B, 
the edge modes become completely flat, apart from a 
small finite-size gap around the zone boundary. We have 
already seen such fiat edge modes in graphene in the 
case of zigzag edge geometry (but on a hexagonal lat- 
tice) . [tI, [111 In graphene, such fiat edge modes connect 
ID projection of K and K' points via the ID BZ bound- 
ary. This is a similar behavior to the present case, if one 
regards the former as the limit of vanishing intrinsic cou- 
pling (or topological mass A) in the KM model. One of 
the differences between the two cases is that here the flat 
edge modes cover the entire ID Brillouin zone. 

In order to elucidate the nature of flat edge mode at 
A — AB, first notice that at this value of A the diagonal 
terms of (the diagonal blocks of the Hamiltonian matrix 
in) Eq. (I43|) vanish. This implies, as in graphene nano- 
ribbon in the zigzag edge geometry, the existence of an 
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FIG. 9: Energy spectrum in the ta^;l&,% edge geometry for different values of A (^4 = B = 1). Upper-left panel: A = 0.25 
(spectrum shown in red), upper-central: A = 0.8-B (spectrum in green), upper-right: A = 2B (spectrum in blue), lower-left: 
A = 3.25 (spectrum in cyan), and lower-central: A = 4B (spectrum in magenta). At A = 4B, the edge modes become 
completely flat and covers the entire Brillouin zone (cf. case of graphene in the zigzag edge geometry 0,112). Notice that here 
the horizontal axe is suppressed to make the edge modes legible. Note that at A = 45 the bulk spectrum is also gapless; the 
completely flat edge modes indeed touch the bulk continuum at the zone boundary (projection of 2D Dirac cones). Compare 
this panel with the corresponding panel of straight edge case: FIG. ((5|. The number of rows Nr is here chosen to be A'^,, — 100. 
These five plots are superposed in the lower-right panel to show that the edge spectra at different values of A are, in contrast 
to the straight edge case, not on the same curve. Even in the long-wave-length limit: — ^ 0, their slopes are different. 



eigenstate of the form: 



be written explicitly as. 





■03 



V'5 



(48) 



i.e., an eigenvector satisfying xp2j — (j = 1,2,---). 
Here, a semi-infinite geometry is implicit (for approxi- 
mating a ribbon of sufficiently large width or AV; our 
system extended from j = 1 to j = Nr), with a bound- 
ary condition of ipo — 0. Under this setup, and with 
the condition of vanishing diagonal matrix elements, Eq. 
P5)) implies, 



rV'sj-hi + rtV'aj-i = (j = i,2,...). 



(49) 



Simultaneously, it should also have a vanishing eigenen- 
ergy i? = for consistency. 

Clearly, Eq. has a solution of the form of a geo- 
metric series, here, for ip2j~i [j — Ij 2, 3, • • • ): 



V': 



2j + l 



(50) 



Note that here A plays, roughly, the role of p^, but their 
precise relation will become clearer when the entire prob- 
lem is solved. In order to proceed, we recall that T can 



r = 



2Bck 
(1 -I- i)(ck + Sk) 



-4(l-i)(cfe-Sfe) 
^2Bck 



(51) 



Then, by assuming a solution of the form of Eq. ()50|) , Eq. 
(|49l) can be reduced to the following eigenvalue problem 
for Tpi: 



with the eigenvalues, 



^" i2cl-l)A^-8B^cl 
and the corresponding eigenvector, u±, i.e., 

— r^"'"r^M± = \±u±, 

given explicitly as, 



a±{l - i) 
1 



(52) 



(53) 



(54) 



(55) 



The coefficient a± is a function oik, which takes precisely 
the following form: 



a±{k) 



-J- f A tan fc T VA^ + 8^2 tan^ k) . (56) 
AB \ / 
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Notice that a_ = — l/(2a+), and the two eigenspinors 
u± are orthogonal: u'_u^ = 0. A general solution in the 
form of Eq. (gS]) can be thus constructed by applying 
— r^^r^, recursively, to 

■01 = + C-U-, (57) 

and the result is, 



a+(l — i) 






1 




1 



(58) 

In this construction, the two eigenvectors u± always have 
the form of Eq. ([55]). This feature remains when A is 
away from 4B at which the edge modes are no longer 
completely flat, or rather even in the regime in which the 
edge spectrum is not flat at all. 

Another remark, concerning the behavior of Eq. ([55)1 
is that under the choice of signs in Eq. ([53)) . |A+| is 
always larger than 1, whereas | A_ | < 1 except at the zone 
boundary. This can be easily verified either numerically, 
i.e., by plotting X± as a function of k, or by showing A— = 
1/ A+ using directly the expression for A± in Eq. (l53l) . In 
numerical simulation for systems of a finite number of 
rows, both of these two solutions play a role giving rise 
to a pair of edge solutions. iSlj 

In FIG. [TOl the upper panel shows A_ plotted as a 
function of fc at A = 1 but for different values of B, 
naturally assuming A = AB. When B is smaller than a 
critical value Be = 0.35... A_ changes its sign (has a zero) 
at intermediate k. [321] This continues to be the case even 
in the limit B vanishes. The lower panel shows the edge 
spectrum at A = 45, i? = 0.2 and A = I. Note that 
the zero of A_ corresponds to the value of k at which the 
bulk spectrum focuses onto a single point. 



1.0 rr 
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•< 0.0 - 
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FIG. 10: A- plotted as a function of at ^ = 1 but for 
different values of B, i.e., B = 1,0.6,0.4,0.35,0.2 and B = 
0, corresponding, respectively, the colors: red, green, blue, 
cyan, magenta and orange (upper panel). Edge spectrum at 
A = 4B and B = 0.2 [A = 1, Nr ^ 100, lower panel). Two 
reference lines are at fc = ±0.438977..., the value of k at which 
A_ vanishes at S = 0.2. The axes are suppressed in the lower 
panel so as to highlight the completely flat edge modes. 



function!) has the following particular form: 



B. Wave functions in special cases of parameters 

As we will describe in detail in the next subsection, our 
"recipe" for constructing the exact edge wave function, 
and simultaneously its spectrum, lies in "extrapolating" 
the exact solution available at A = 4i? to a general value 
of A/B € [0,8] (recall also FIG. (J). To complete this 
program, we need to refer to some results of the numerical 
experiments performed for a system of finite number of 
rows. We have already seen the spectrum of such systems 
in FIG.ini here we focus on the behavior of wave function, 
i.e., the behavior of ipj as a function of j. 

In the zigzag edge geometry, it is remarkable that (not 
only) the edge wave function (but also the bulk wave 



ci(l-i) 

C2 

C5(l - l) 
C6 



(59) 



when that eigenstate represents an edge mode, Eq. 
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further simplifies: 



C2 



C4 



C6 



a{l-i) 
1 

a{l-i) 
1 

a{l~i) 
1 



(60) 



i.e., for a given set of parameters A, B (and A) as well 
as for a fixed k, Uj = C2j-i/c2j is a constant {— a). The 
ratio, on the other hand, 



Pj 



C2j+2 
C2j 



(61) 



is a measure of, to what extent the edge mode is localized 
in the vicinity of a boundary, say, at j = 1. 




FIG. 11: The ratio, pj = C2j+2/c2j {j = 1,2, 3, •■•) plotted 
(red points with filling to the horizontal axis) at A: = 0.05 for 
A/B = 0.25 (upper panel) and A/_B = 0.30 (lower panel). 
A = B = 1, Nr = 100. In the upper panel, the blue 
line corresponds to the "theoretical" value, pj = 0.759908..., 
whereas in lower panel, the plots are fitted by a curve, 
pj — rsin(j + 1)9/ sin jd, with the choice of parameters, 
r = 0.691189..., 9 = 0.133449... 




FIG. 12: Theoretical value (derived later) of p (its magnitude, 
\p\) plotted as a function k for A/B = 0.25 [A = B ^ 1). At 
kf-K — 0.05 there are two possible solutions for p: p_i = 
0.759908... and p_2 = 0.628683..., the larger value of which 
determines large-j behavior of pj — C2j+2/c2j. Merger with 
bulk occurs when p_i — 1, i.e., at kf-K = 0.263808... Close 
to the zone boundary {k/n > 0.900237...), reentrance of edge 
solution occurs (see FIG. [17] for details). 




FIG. 13: Same as FIG. [H for A/B = 0.3 {A = B = 1). 
On the reference fine at = 0.05, p ~ 0.685038 ± 0.0919993i 
(IpI = 0.691189..., Arg[p] = ±0.133499...). 



We have extensively studied such characteristic behav- 
iors of the edge wave function in numerical experiments. 
In FIG. [Til results of such analyses are shown, for the 
choice of parameters such that k/ir = 0.05, A = B = 1 
and two different values of A/B: A = 0.25 and A = 0.30 
for comparison. At this value of k/n — 0.05, we first ver- 
ified that the wave function ipj takes indeed the form of 
Eq. (|60p. with a± approximately given by. 



a+ = 0.687705..., a_ = -0.727056. 



(62) 



As in the straight edge case, for k corresponding to an 
edge mode, i.e., for such a state that are localized in 
the vicinity of either of the two boundaries, this value 
of a± is common practically to all j in the strip, as far 
as a finite amplitude exists. For bulk states which are 
well extended into the interior of the sample, the wave 
function takes no longer the form of Eq. (|60p , but keeps 
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still a characteristic form as Eq. (1591) . 

In the two panels of FIG. [TTJ pj is plotted as a function 
of j for two different values of A/B. At A = 0.25, pj 
saturates at rows away enough from the boundary at j ~ 
1 (but not too close to the other edge, either), to a value 
close to Pj = 0.759908... (blue line), a value which is later 
"derived" (see FIG. [H]). Let us assume, as in the 
straight edge case, that the wave function in the zigzag 
edge geometry takes the following form: 



where the eigenspinors u± are always given by Eq. (1551) , 
the latter found analytically in the limit A/B 4. We 
have extensively verified the validity of this hypothesis 
in numerical experiments. The coefficients c±i^2 are sus- 
ceptible of system's geometry. Here, in a strip geometry, 
they satisfy. 



c+i = 7^ 0, C_i = C_2 

for one edge mode, and 



0. 



c+i = C+2 = 0, c_i = -C_2 ^ 



(63) 



(64) 



for the other. Under this hypothesis, such a behavior 
as seen in the upper panel of FIG. [11] (A — 0.25) is 
interpreted as a consequence of two "real solutions" for 
p, which are also both smaller than unity (cf. FIG. [T^ . 

On the other hand, at A = 0.30 (in the lower panel 
of FIG. [TT]) Pj shows an oscillatory behavior. This im- 
plies, with the same hypothesis as above, i.e., the wave 
function, ij^ji given as in Eq. (|63l) . with the choice of 
coefficients as Eqs. (j64p . a pair of complex solutions for 
p. Indeed, the plots at A = 0.30 are nicely fitted by a 
curve, of the form, 



P3 



| Sin(j + l)6' 
sin j0 



r [cos6l-hsin6'cot(J0)] , (65) 



with the choice of parameters, \p\ = 0.691189... and = 
0.133499..., which will be also (a posteriori) justified (see 
FIG. El). 

Comparing these two contrasting cases, notice that the 
coefficients a±, estimated to be such as Eq. (p^ . are 
common to the two cases, i.e., independent of A/B. One 
can indeed verify (by changing the parameter A/B in 
numerical experiments) that the eigenspinors u± remains 
always the same, as far as the state describes an edge 
mode (see FIG. [14]); only pj changes as a function of 
A/B. 

In FIG. [TH aj = C2j-i/c2j (j = 1, 2, 3, • • • ) is plotted 
for the lowest-energy eigenmode ^I^o (in the upper band) 
at different values of A/B. One can see that in the range 
of k at which 'J'o is expected to represent an edge mode 
the plotted points fall roughly on the theoretical curve for 
a_ (fc) — cf. Eq. (j56l) — apart from a small disagreement 
close to the zone boundary (fc = tt). This is indeed a 
key discovery allowing us to proceed to the next step, of 
extrapolating the earlier exact solution at A = AB to an 
arbitrary value of A/B. 




FIG. 14: aj = C2j-i/c2j {j = 1,2,3, •■•) calculated in a 
strip geometry (here, Nr = 100) is plotted against the the- 
oretical curve for q_ (k) (shown in black) given in Eq. (|56p 
up to j = Nr/4: for the lowest-energy eigenmode at differ- 
ent values of A/B = 0.1, 0.25, 0.5, 1, 1.35 (A = B = 1 fixed), 
corresponding, respectively, to the colors: red, green, blue, 
cyan and magenta. Reference lines indicate the regime of k in 
which the edge modes disappear at the given value of A/B. 



C. Derivation of exact edge wave functions 

Let us reformulate the recipe for constructing the exact 
edge wave function and simultaneously its spectrum in 
the zigzag edge geometry, which has already been briefly 
outlined in the introduction (recall also FIG.[T|). 

1. We have seen in the previous subsection that the 
edge wave function \1/ in the zigzag edge geometry 
always takes, as far as it describes a localized edge 
mode, the form of Eq. [60] with a parameter a± 
depending only on k (and A, B). All our numer- 
ical data agree with the hypothesis that a±, con- 
sequently the reduced two-component eigenvector 
u±, is independent of A. 

2. On the other hand, we know that the problem can 
be solved exactly at A = AB. We have seen, in 
particular, that the wave function ^I^ can be con- 
structed from the same set of spinors u± with a 
choice of parameters a± given analytically as a 
function of k in Eq. ([56]) . 



Taking also into account the fact that the edge modes, 
gapless at fc = and characterizing the topological in- 
sulator, evolves continuously to the completely flat edge 
mode at A = 4i?, one can deduce, from these two obser- 
vations, that the solution of the eigenvalue equation for 
"00 J he., Eq. (|45p for an arbitrary A should be given, in- 
deed, by u±, deflned as in Eq. (|55p . with the parameter 
a±(fc) obtained analytically in the limit: A = AB (recall 
FIG.[14]). 

Thus, for a general value of A, only p and E are un- 
known (recall that the edge spectrum is no longer flat for 
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kin kin 

FIG. 15: Four solutions for p (each curve corresponds to one solution). Only the magnitude of such solution, which is generally 
a complex number, i.e., |p| is plotted as a function k at different values of A/B = 0.25, 0.3, 1.35 and 1.45 {A = B = 1). 



irt' 




a{\ - i) 


= E 


a{l - i) 


p . 




1 




1 



(66) 



a general A). But, clearly, they are solutions of 



(A - 4B)ct^ + pT 



where the 2x2 matrix T is given explicitly as Eq. ([5T|). 
34| To find exactly the value of p and one has only 
to solve this set of equations, and at the end of the cal- 
culation, substitute the analytic expression for a, i.e., 
a± given in Eq. ((56)) obtained in the limit of A = AB. 
Clearly, Eqs. (|66|) are a set of coupled equations, linear 
in E and quadratic in p. We expect, therefore, two sets 
of solutions for (p, E) , which are given as a function of 
a. To each of these two sets of solutions, one substitutes 
either a = a+ or a = There exist, therefore, four 
sets of solutions, in general. 

Unfortunately, the analytic formula for these four sets 
of solutions are too lengthy to be shown here. Instead, 
we plotted these four solutions in FIG. [TH for different 
values of A/B. 



D. Reentrant edge modes 

Reentrance of the edge mode is another characteristic 
feature of the edge mode of zigzag geometry, and oc- 
curs close to the zone boundary, k/n = 1, when A/B 
is not too large: A/B < 1.354.... Very remarkably, the 



spectrum looks completely "innocent" when this occurs, 
i.e., the edge mode, say, the lowest energy {— Eo) mode 
in the upper band looks almost completely degenerate 
with the bottom of the (bulk) spectrum {— Ei), in this 
regime of k (see FIG. [18]). Existence of an edge mode 
of such specific character is, on the other hand, noth- 
ing exceptional in the zigzag edge geometry. At a value 
of A, e.g., A/B = 0.25 or A/B = 0.30 as in FIG. HH 
such reentrant edge modes are indeed existent. If one 
focuses on the wave function of, say, the lowest energy 
mode in the upper band, after touching the lower band 
at fc = 0, it continues to be spatially localized when k 
is small enough, but as the spectrum merges with the 
bulk continuum, the wave function also starts to pene- 
trate into the bulk. However, close to the zone boundary, 
it starts to be localized again. This is what we call the 
reentrance of edge modes. 

Figs, [in] and [T7] highlight the behavior of such reen- 
trant edge modes, naturally in a different regime of k 
from, say, FIG. [11] Fig. [1^] shows the behavior of pj 
at k/n = 0.901 and fc/vr = 0.903 for A = 0.25. At this 
value of k, the wave function ipj takes always the form 
of Eq. (|5(I)) . but with a different set of parameters for 
a± from the case of FIG. [TT] upper panel, since it still 
depends on k. The two plots for pj in Fig. [TBI show two 
typical behaviors of the edge wave function, i.e., one cor- 
responding to real and the other to complex solutions for 
p. As we have extensively studied in the case of ordinary 
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FIG. 16: Reentrant edge modes I: pj plotted as a function 
of j = 1, 2, 3, • • • at A = B = 1 and A = 0.25 (as in FIG. 
111! upper panel), but for k/-n = 0.901 (upper panel) and 
k/n = 0.903 (lower panel). In the upper panel, the blue line 
corresponds to the theoretical value, p — 0.956432 whereas, 
in the lower panel the plots (in red) are fitted by the curve, 
Pj = rsin(j + 1)6/ sin j9, with the choice of parameters, r = 
0.867804, e = 0.140687. 
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FIG. 17: (Theoretical value of) \p\ plotted as a function k 
at A = 0.25 (again, as in FIG. [TTl upper panel). Enlarged 
picture for k close to the zone boundary: fc/yr ~ 0.9 The reen- 
trance of edge solution occurs at k/n = 0.900237..., whereas 
two real solutions for p is possible when k/n < 0.901675... 
Two reference lines are at fc/vr = 0.901 and a.t k/n — 0.903, 
on which pj was plotted in FIG. 1161 The value of p on these 
lines are, p = 0.783429... and p = 0.956432... on k/n = 0.901 
(case of real solutions), whereas p ~ 0.859242 ± 0.121601j on 
k/n = 0.903. 




FIG. 18: Merger of the edge mode with bulk continuum and 
"absence" of reentrance in the spectrum. As for the latter, 
it turns out later that binding energy of the edge mode is 
too small to be seen at this scale (see FIG. I19p . This is an 
enlarged image of the spectrum in the zigzag edge geometry 
as shown in FIG. (9] Here, the parameters are chosen such 
that A/B = 1.2, A = B ^ 1, Nr ^ 100. Two reference lines 
at k/n = 0.643658... = kci/n and at k/n = 0.826568... = 
kdi/n introduce three different momentum regions: i) k/n £ 
[0,kci/n], ii) k/n £ [kci/n,kc2/n] and iii) k/n e [kc2/n,l], 
corresponding, respectively, to i) the ordinary edge, ii) the 
bulk and iii) the reentrant regimes. 



edge modes (appearing at /c/tt <C 1), the two contrasting 
behaviors of pj as a function j (in FIG. fTTj) are naturally 
understood by referring to the theoretical curve of p as a 
function of k, e.g., such as the one shown in FIG. [T^ 

What is rather remarkable here, in the case of Fig. [T51 
is that this crossover between real and complex solutions 
occurs within a tiny change of k, i.e., from k/ir = 0.901 
in the upper panel to k/n = 0.903 in the lower panel. 
This drastic change is, however, quite reasonable from 
the viewpoint of FIG. [T71 The upper panel of Fig. [12] 
shows a monotonic decay, which converges asymptoti- 
cally to a single exponential decay. This is consistent 
with the behavior of theoretical curve for p as a function 
of k in FIG.[T71 The latter implies two real solutions for 
p at k/n = 0.901: p = 0.783429... and p ^ 0.956432... 
The latter coincides with the value of pj in FIG. [12] 
at which it saturates. At k/n — 0.903, on the other 
hand, the plots for pj are nicely fitted by the curve, 
Pj = rsin(j + 1)0/ sin with the choice of parameters, 
r = 0.867804... and 9 = 0.140687... (see the lower panel 
of Fig. [TO)) . This is a clear fingerprint that the reentrant 
edge mode at this value of k corresponds to a pair of 
complex solutions for p. 

Does the reentrant edge mode really have zero binding 
energy? In order to address this question, we (re)plottcd 
the energy spectrum {Ei — Eq, to be precise) but in an 
enlarged scale roughly by one thousand times in FIG. 
[T9l First, for an "ordinary" edge state, occurring at 
< \k\/n < 0.289936 = k^/n the value Ei - Eq is 
much above the threshold at this scale. Ei — Eq takes a 
value of order ~ 1 for such ordinary edge state. As for 
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a 0.004 




FIG. 19: Binding energy of the reentrant edge mode: Ei ~ Eo 
is plotted as a function of fe at A = 0.3B, A = B — 1. 
Different curves correspond to different size (width) of the 
system: iV^ = 100 (blue), 200 (green) and 300 (red). The 
two reference lines are placed at fc/vr = 0.289936 = kci/n and 
k/-K = 0.8983 = fcc2/7!" (fc G [kci,kc2] corresponds to the bulk 
regime). The plots reveal an extremely small but still a finite 
binding energy of the reentrant edge modes. 



the reentrant edge mode, FIG. [12] reveals that it has in- 
deed an extremely small but still a finite binding energy. 
Notice different behaviors of Ei — £^o as a function of Nr 
in the bulk and reentrant regions of k. The former (the 
latter) corresponds to k/n e [kci/TT,kc2/'^ = 0.8983...] 
{k/n € [A;c2/7r, 1]). In the bulk region Ei — Eq is ex- 
pected to vanish in the thermodynamic limit. FIG. \W\ 
shows indeed that the binding energy of reentrant edge 
mode. El — Eq ^ 0.001, is thousand times smaller than 
that of the ordinary edge state. This implies the appear- 
ance of an extremely small energy scale which was not 
existing in the original Hamiltonian (cf. Kondo effect). 

The reentrance of edge mode is indeed a unique feature, 
in its contrasting properties in real and momentum space, 
and in the appearance of an extremely small energy scale. 



VI. CONCLUSIONS 

We have highlighted in this paper various unique prop- 
erties of helical edge modes in Z2 topological insulator. 
We have extensively investigated a lattice version of the 
BHZ model, under different edge geometries. One of the 
specific characters of BHZ model is that the spin Hall 
conductance in the bulk changes its sign in the middle 
of topological phase (at A = 4i?), i.e., aiy = ±e/(27r), 
respectively, for < A/B < 4 and for 4 < A/B < 8, 
though both represent a non-trivial value. From the 
viewpoint of bulk-edge correspondence, this information 
should be also encoded in the edge theory. We have seen 
that the change of ax^J manifests in a very different way 
in the (1,0)- (straight) and (1,1)- (zigzag) edge geome- 
tries. In the (l,0)-edge case, the edge spectrum changes 
its global structure in the two parameter regimes, i.e., 
the main location of the mode moves from the zone cen- 



ter for A/B < 4, to the zone boundary for A/B > 4. 
As a result, the group velocity at the intersection with 
Fermi level reverses its sign, leading to change of the sign 
in Landauer conductance at the edge. In the (1, l)-edge 
case, on the other hand, the edge spectrum is symmetric 
w.r.t. A — AB, i.e., neither change of the position of gap 
closing, nor the reversal of group velocity at A = AB. 

(s) . 

The change of axy is here encoded in the swapping of 
left- and right- going edge modes of the same spin. 

Much of our focuses has been on the analysis of the 
zigzag or (l,l)-edge geometry, the latter showing, as a 
consequence of specific way in which the bulk topologi- 
cal structure is projected onto the ID edge, a number of 
unique features, such as the completely fiat edge spec- 
trum at A = 4B, and the reentrance of edge modes. We 
have also shown, here in a half-empirical way, that the ex- 
act edge wave function for zigzag edge geometry can be 
constructed, by extrapolating the solution at A = AB. 
The reentrant edge mode, though sharing much of its 
characteristics with the usual edge mode in real space, 
introduces a new extremely small energy scale which was 
absent in the original BHZ model. 
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Appendix A: Edge solution in the long-wave-length 
limit [12, 24, 25] 

Let us first recall that the eigenvector |d(fc)±), which 
has appeared in Eq. (jS]), is a standard SU(2) spinor, here 
chosen to be single-valued. An eigenvector, correspond- 
ing to a positive energy eigenvalue > 0, is \d{k)+) 
with k satisfying E — E{k), and represented as. 



|d(fc)- 



e-"^ cos(6'/2) 
sin(6l/2) 

1 

^/2did-d,) 



(Al) 



dx idy 
d~ d. 



9, (j) are polar coordinates in d-space, satisfying the rela- 
tions such as, 



cos9 — coscj)— — ^ (A2) 
d Idl + dl 
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In order to find an edge solution, we focus on a solution 
of the form, pJl] 



(A3) 



say, in a semi-infinite plane: y < 0. The spatial depen- 
dence in the y-direction can be taken into account by 
applying Pier Is substitution: ky — >■ —id/dy to fcy's in 
h{k). 

The eigenenergy of such a solution is obtained by a 



simple replacement: ky 



in in Eq. ([8|), i.e.. 



= A^ + {A^ - 2BA) {kl - n^) + B^{kl ~ n^f. (A4) 

This can be regarded as a quadratic equation w.r.t. k^. 
Its two solutions are. 



- 2BA 
2^2 



± — [A^ - 4BA) + 4B^E^. (A5) 

For a given set of kx and there are two possible val- 
ues for K^, or equivalently, four possible values for k. Of 
course, in a semi-infinite plane, say, y < the edge solu- 
tion of the form, Eq. (IA3|) should decay as ?/ — > — cx), so 
only two of such solutions are relevant. 

We expect that the edge spectrum behaves as in 
the limit offc^; — >■ 0. Let us parametrize the two solutions 
in this limit as 



,(0) 



P±i 



(A6) 



where 



P 



A^ - 2BA 
2^2 ' 



Q^—^A^ {A^-4BA). (A7) 



When < A < A'^/{4B) = Ai, P > and Q is real. 
Since \a\ > 6 as far as 5 is real, Eq. (jASP represents two 
positive solutions for k^, i.e., the wave function repre- 
sented by Eq. (lASp shows simple exponential damping. 
On contrary, we expect a pure imaginary solution for k 
for an extended state in the bulk. On the other hand, 
when Ai < A < Aq = A'^/{2B), P is always positive 
but Q becomes purely imaginary. Thus two solutions 
for K become complex numbers conjugate to each other. 
In this case, the wave function represented by Eq. (IA3[) 
shows damped oscillation. 

The corresponding eigenvector is obtained by the same 
replacement ky —in, here in Eq. (jA2p . i.e.. 







A(kx ^) 






_E-d,{K) _ 



\d{kx,-iK)- 



(A8) 

where ^^(k) = A + B{k^ — k^). For a given value of k^ 
and -E > 0, we thus have identified two solutions charac- 
terized by different values of k. In order to construct a 



general solution in the presence of a boundary, we need 
to take a linear combination of these two solutions, i.e.. 



(A9) 



where u±, v± are short-hand notations for and 

We now fix the boundary condition at y = 0, which we 
choose to be. 



(AlO) 







c+ 











c_ 








which implies the following secular equation 



det 



A{kx - K+) A{kx - K-) 

E-d,{K+) E-d,{K-) 



= 



This leads to, 

E{kx) = A - Bk+k^ + Bkx{n+ + k_) - Bkl. 



(All) 



(A12) 



We have thus identified the two basic equations, Eqs. 
(|A5[) and (jA12[) . for determining the energy spectrum 
E = E{kx). 

Let us check whether this solution contains the edge 
modes. We expect that the edge spectrum bahaves, as 
kx ^ 0, E ^ 0. Recall that in this limit, Eq. (|X5|) 
reduces to Eqs. (|A6I) and (IA7|) . Eq. ([57)1 is also simplified 
in this limit, as 



E = E{0)^A-Bk^"^k^°\ 



(A13) 



Focusing on the case, k±(0) > and B,A > 0, and 
using the parameterization in Eqs. (|A6p and (|A7p . one 
can readily verify. 



(0) (0) 



(A14) 



Thus Eq. (IA13|) is safely satisfied. 

How about the first order corrections? i.e., contribu- 
tions of order 0{kx) to the energy spectrum, E = E{kx)- 
First note that there is no C'(A:a;)-correction to k±. One 
can, therefore, safely replace, at this order, k±'s in Eq. 
([57)1 with their values at fc^; — > 0, i? — ?> 0, i.e.. 



E = E{0) = A - 



Bkx 



AO) , ^(0) 



(A15) 



We have already seen that the first two terms cancel, 
whereas 



,(0) 



,(0) 



2a + 2VP2-Q2 



(A16) 



Thus, the edge spectrum in the continuum limit is deter- 
mined to be, 



E{kx) = ±Akx + 0{kl 



(A17) 



19 



Remarkably, the slope of the edge spectrum depends only 
on a single parameter, A. An interesting question is to 
what extent this conclusion is general? If one calculates 
the edge spectrum, using a tight-binding model, generally 
the results depend on the way edges of the sample are 



introduced with respect to the lattice. In the case of 
zigzag edge, in particular, apparently the edge spectrum 
does not converge to Eq. ()A17P even in the long-wave- 
length limit: fcj, — ^ 0. 
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integer multiple of 2e^ /h in the language of charge con- 
ductance. 

In the case of graphene (and also KM) zigzag edges, we 
adopt a different boundary condition: only the A (B)- 
sublattice component of the wave function vanishes at 
one (the other) boundary. For details, see, Ref. 0, [3]- 
The bulk solutions of Eq. (|22|) corresponds to the choice, 
p — g'^^y^ or IpI = 1 which is consistent with the Bloch 
theorem. In the strip geometry with the periodic bound- 
ary condition, ky takes discrete values. With the open 
boundary condition relevant to the actual strip, ky is no 
longer a good quantum number. However, if the width 
of the strip is large, one may roughly interpret the ID 
energy spectra in the strip geometry as composed of the 
many slices of bulk energy spectrum at different values 
of ky. In addition, a pair of edge modes appear as a char- 
acteristic feature of the nontrivial topological property. 
Real solutions for p appear in the regime: < A < A_, 
and also at the other end. In the approximation that 
becomes valid in the small wave number, this threshold 
value is given by Ai = A^/(4B). For A = B = 1, A_ = 
2-73 = 0.2679 • • • , whereas, Ai is, of course, 1/4 = 0.25. 
If one expands Eq. (|4ip in powers of ^^/(4_B^), then at 
leading order A_ coincides with Ai. 
Namely, if one compares it to the straight edge case, 
e.g., in considering the long-wave-length limit, one has 
to make the correspondence between ky/2a and k^a. 
When Nr is odd, one solution, corresponding, say, to 
A_, is localized in the vicinity of ji = 1, and indeed has 
the form of Eq. H48|) with '(/'2i+i given as Eq. (|58p and 
c+ — 0. The other solution, corresponding to A+, has the 
same structure of Eq. (|48|) but with tp2j+i increasing with 
practically an equal geometric ratio of A+, and naturally 
localized in the vicinity of the other edge: j = Nr. On 
the other hand, when A''^ is even, the eigenmode of the 
system becomes a linear combination of the above two 
types of solutions. This even/odd feature occurs only at 
A precisely equal to 4_B, since at this value of A where 
the edge spectrum becomes completely flat, the two (gen- 
erally) counter-propagating edge modes acquire the same 
(zero) group velocity, and get mixed. 
Clearly, at this value of k the edge wave function is ex- 
tremely localized, i.e., onto a single row: j = 1 or j = Nr. 
We leave formal derivation of Eq. (|63|) to a future publi- 
cation. Instead, we take it here as a plausible hypothesis 
fully justified by numerical experiments. 
Inspecting the explicit form of Eq. (|66|) and (|51|l , notice 
that 1 ±i factors out. So all the coefficients become real. 
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